################################################################################
#############Authoritarian footprints in Central and Eastern Europe#############
###############################Replication Script###############################
################################################################################

##########LIBRARIES AND FUNCTIONS##########
library(plyr)
library(dplyr)
library(zoo)
library(broom)
library(lubridate)
library(grid)
library(gridExtra)
library(foreign)
library(stargazer)
library(lattice)
library(lme4)
library(gsubfn)
library(reshape)
library(ggplot2)
library(RColorBrewer)


###extracting legend of plot###
g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}#function to extract legend from plot
###creating shared legend for multiple plots###
grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {
  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + 
                    theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x +
                 theme(legend.position = "none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)
  
  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl), 
                                            legend,ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend, ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))
  
  grid.newpage()
  grid.draw(combined)
  
  # return gtable invisibly
  invisible(combined)
}#function to arrange plots
###country cluster test formula
get_confint<-function(model, vcovCL){
  t<-qt(.975, model$df.residual)
  ct<-coeftest(model, vcovCL)
  est<-cbind(ct[,1], ct[,1]-t*ct[,2], ct[,1]+t*ct[,2])
  colnames(est)<-c("Estimate","LowerCI","UpperCI")
  return(est)
}
cluster.se<-function(model, cluster)
{
  require(multiwayvcov)
  require(lmtest)
  vcovCL<-cluster.vcov(model, cluster)
  coef<-coeftest(model, vcovCL)
  return(coef)
}



##########DATA IMPORT##########
setwd("DIRECTORY WHERE DATA IS SAVED / GRAPHS ARE EXPORTED TO")
af_main <- read.csv("./authoritarian_footprints_replication_main.csv")
af_wide <- read.csv("./authoritarian_footprints_replication_wide.csv")


##########FIGURE 1##########

###defining outcome functions
functions <- c("INDLIB", "RULEOFLAW", "PUBLIC", "COMPET", "MUTUCONS", "GOVCAP", "TRANSPAR", "PARTICIP", "REPRES")
functions_var <- c("indlib","ruleoflaw","public","compet","mutucons","govcap","transpar","particip","repres")

###get data into right format
rtrend <- af_main[c("country","year","region2",functions)]
for (i in (functions)) {
  step1 <- paste("rtrend <- rtrend %>% group_by(country) %>% arrange(country,year) %>% mutate(",i," = na.approx(",i,", na.rm=FALSE, rule =2))",sep="")
  step2 <- paste("rtrend <- rtrend %>% group_by(country) %>% arrange(country,year) %>% mutate(",i," = na.locf(",i,", na.rm=FALSE))",sep="")
  step3 <- paste("rtrend <- rtrend %>% group_by(country) %>% arrange(country,year) %>% mutate(",i," = na.locf(",i,", na.rm=FALSE, fromLast=T))",sep="")
  eval(parse(text=c(step1,step2,step3)))
}
for (i in (functions)) {
  step1 <- paste("rtrend <- rtrend %>% group_by(region2,year) %>% arrange(region2,year) %>% mutate(",i," = mean(",i,",na.rm=T))",sep="")
  eval(parse(text=c(step1)))
}
rtrend$region2 <- as.factor(gsub(" ", "_",rtrend$region2))
rtrend <- data.frame(unique(rtrend[c("year","region2",functions)]))
rtrend <- melt(rtrend, id=c("region2","year"))
rtrend$variable <- revalue(rtrend$variable, c("INDLIB"="Individual Liberties", "RULEOFLAW"="Rule of Law", "PUBLIC" = "Public Sphere", "COMPET" = "Competition", "MUTUCONS" = "Mutual Constraints", "GOVCAP" = "Governmental Capacity", "TRANSPAR" = "Transparency", "PARTICIP" = "Participation", "REPRES" = "Representation"))

###plot
colours <- brewer.pal(9, "Set1")
colours[6] <- "#000000"
figure1 <- ggplot() + 
  geom_line(data = rtrend, aes(x=year, y=as.numeric(paste(value)), colour=variable, linetype=variable)) + 
  scale_color_manual(values = colours) + theme_bw() + xlim(1989, 2017) + ylim (0, 100) + 
  labs(title = 'Central_and_Eastern_Europe', x ='Year', y = 'Standardized Value', colour = 'Democratic Principle', linetype= 'Democratic Principle') + 
  theme(text=element_text(family='Times New Roman')) + 
  ggtitle("Central and Eastern Europe")
ggsave(file="./figure1.png", figure1, width = 8.27, height = 3, dpi = 600)


##########FIGURE 2##########

###get data in right format I: country trends
ctrend <- af_main[c("country","year",functions)]
ctrend$country <- gsub(" ", "_",ctrend$country)
ctrend$country <- as.factor(gsub("-", "_",ctrend$country))
ctrend <- subset(ctrend,country !="Kosovo")
levels(ctrend$country)[levels(ctrend$country)=="Serbia_/_Yugoslavia"] <- "Serbia"
countries <- list(levels(as.factor(paste(ctrend$country))))
ctrend <- melt(ctrend, id=c("country","year"))
ctrend$variable <- revalue(ctrend$variable, c("INDLIB"="Individual Liberties", "RULEOFLAW"="Rule of Law", "PUBLIC" = "Public Sphere", "COMPET" = "Competition", "MUTUCONS" = "Mutual Constraints", "GOVCAP" = "Governmental Capacity", "TRANSPAR" = "Transparency", "PARTICIP" = "Participation", "REPRES" = "Representation"))

###get data in right format II: populism periods
ctrend_pop <- af_main[c("country","year","populism_gov")]
ctrend_pop$country <- gsub(" ", "_",ctrend_pop$country)
ctrend_pop$country <- as.factor(gsub("-", "_",ctrend_pop$country))
levels(ctrend_pop$country)[levels(ctrend_pop$country)=="Serbia_/_Yugoslavia"] <- "Serbia"
ctrend_pop$populism_gov <- as.factor(paste(ctrend_pop$populism_gov))
ctrend_pop_l1 <- ctrend_pop
ctrend_pop_l1$year <- ctrend_pop$year + 1
ctrend_pop_l1$populism_gov_l1 <- ctrend_pop_l1$populism_gov
ctrend_pop_l1 <- ctrend_pop_l1[c("country","year","populism_gov_l1")]
ctrend_pop <- left_join(ctrend_pop,ctrend_pop_l1,by=c("country","year"))
ctrend_pop$period <- ifelse(ctrend_pop$populism_gov != ctrend_pop$populism_gov_l1, 1, 0)
ctrend_pop$period <- ifelse(is.na(ctrend_pop$period), 1, ctrend_pop$period)
ctrend_pop <- ctrend_pop %>% group_by(country) %>% arrange(country,year) %>% mutate(period = cumsum(period))
ctrend_pop <- ctrend_pop %>% group_by(country, period) %>% arrange(country,year) %>% mutate(from = min(year))
ctrend_pop <- ctrend_pop %>% group_by(country, period) %>% arrange(country,year) %>% mutate(to = max(year))
ctrend_pop <- unique(ctrend_pop[c("country","from","to","populism_gov")])
ctrend_pop_pop <- unique(ctrend_pop[c("country")])#adds a year 0 to year 0 row, so that the rectangle tool below does not get a NULL dataframe
ctrend_pop_pop$from <- 0
ctrend_pop_pop$to <- 0
ctrend_pop_pop$populism_gov <- "1"
ctrend_pop_na <- unique(ctrend_pop[c("country","populism_gov")])#adds a year 0 to year 0 row, so that the rectangle tool below does not get a NULL dataframe
ctrend_pop_na$from <- 0
ctrend_pop_na$to <- 0
ctrend_pop_na$populism_gov <- "NA"
ctrend_pop <- rbind(ctrend_pop, ctrend_pop_pop, ctrend_pop_na)

###creating individual country graphs
countries <- unlist(as.list(levels(as.factor(ctrend$country))))
for (i in (countries)) {
  step1 <- paste(i," <- subset(ctrend, country =='",i,"')",sep="")
  step2 <- paste(i,"_pop <- subset(ctrend_pop, country =='",i,"' & populism_gov == '1')",sep="")
  step3 <- paste(i,"_pop_na <- subset(ctrend_pop, country =='",i,"' & populism_gov == 'NA')",sep="")
  step4 <- paste(i,"_plot <- ggplot() +
                 geom_rect(data=",i,"_pop, aes(xmin = from - 0.5, xmax= to + 0.5, ymin=-Inf, ymax=+Inf, fill='pink'), colour=NA, alpha=0.2) +
                 geom_rect(data=",i,"_pop_na, aes(xmin = from - 0.5, xmax= to + 0.5, ymin=-Inf, ymax=+Inf, fill='grey'), colour=NA, alpha=0.2) +
                 scale_fill_manual('Populist party in government coalition',values = c('grey', 'pink'), labels = c('No information','Yes'),  
                 guide = guide_legend(override.aes = list(alpha = 0.2))) +                 
                 geom_line(data = ",i,", aes(x=year, y=as.numeric(paste(value)), colour=variable, linetype=variable),size=1) + scale_color_manual(values = colours, guide = guide_legend(nrow = 3)) +
                 theme_bw() + scale_x_continuous(limits = c(1989,2017), breaks=c(1990,2000,2010)) +
                 ylim (0, 100) +
                 labs(title = '",i,"', x ='Year', y = 'St. Value', colour = 'Democratic Principle', linetype = 'Democratic Principle') +
                 theme(legend.box = 'vertical', text=element_text(family='Times New Roman', size = 16))",sep="")
  eval(parse(text=c(step1, step2, step3, step4)))
}

###combining individual graphs, relabeling, and exporting
figure2 <- grid_arrange_shared_legend(Albania_plot,Bosnia_Herzegovina_plot + ylab("") + ggtitle("Bosnia"),Bulgaria_plot + ylab(""),Croatia_plot,Czech_Republic_plot+ ylab("") + ggtitle("Czech Republic"),Estonia_plot + ylab(""),Hungary_plot,Latvia_plot  + ylab(""),Lithuania_plot + ylab(""),Macedonia_plot,Moldova_plot + ylab(""), Montenegro_plot + ylab(""), Poland_plot,Romania_plot + ylab(""), Serbia_plot + ylab(""),Slovakia_plot,Slovenia_plot + ylab(""),Ukraine_plot + ylab(""), ncol=3,nrow=6)
ggsave(file="./figure2.png", figure2, width = 8.27, height = 11.69, dpi = 600)


##########FIGURE 3a / B1a: POPULISTS IN OPPOSITION##########

for (i in (functions_var)) {
  step1 <- paste("popopp_",i," <- subset(af_wide, populism_opp == 1 & !is.na(",i,"))",sep="")
  step2 <- paste("popopp_",i," <- popopp_",i,"[c('gwid','gwid_populism_opp_y','populism_opp_y','",i,"')]",sep="")
  step3 <- paste("popopp_",i,"$",i,"_gm <- mean(popopp_",i,"$",i,", na.rm=T)",sep="")
  step4 <- paste("popopp_",i,"$populism_y_gm <- mean(popopp_",i,"$populism_opp_y, na.rm=T)",sep="")
  step5 <- paste("popopp_",i," <- popopp_",i," %>% group_by(gwid_populism_opp_y) %>% mutate(",i,"_periodm = mean(",i,",na.rm=T))",sep="")
  step6 <- paste("popopp_",i," <- popopp_",i," %>% group_by(gwid_populism_opp_y) %>% mutate(populism_y_periodm = mean(populism_opp_y,na.rm=T))",sep="")
  step7 <- paste("popopp_",i,"$",i,"_demean <- popopp_",i,"$",i," - popopp_",i,"$",i,"_periodm + popopp_",i,"$",i,"_gm",sep="")
  step8 <- paste("popopp_",i,"$populism_y_demean <- popopp_",i,"$populism_opp_y - popopp_",i,"$populism_y_periodm + popopp_",i,"$populism_y_gm",sep="")
  count <- paste("popopp_",i,"$p10 <- ifelse(nrow(unique(popopp_",i,"[c('gwid')])) >= 10, 1, 0)", sep="")
  step9 <- paste("popopp_",i,"_plot <- ggplot() +
                 geom_line(data = popopp_",i,", aes(x = populism_opp_y, y = ",i,", group = gwid_populism_opp_y)) +
                 stat_smooth(data = popopp_",i,", aes(x = populism_y_demean, y = ",i,"_demean), method = 'lm', fullrange=TRUE, xseq = seq(0, max(popopp_",i,"$p10 * popopp_",i,"$populism_opp_y, na.rm=T), length=80), col = 'red') + 
                 labs(x ='Populism in opp. (y.)') +
                 theme(legend.position='none',plot.margin=unit(c(0.25,0.5,0.25,0.25), 'cm')) +
                 coord_cartesian(ylim = c(0, 100), xlim =c(1,10)) + 
                 scale_x_continuous(expand = c(0,0), lim = c(0,10), breaks = c(2,4,6,8,10)) +
                 theme(text = element_text(size=rel(3.25))) +
                 theme(text=element_text(family='Times New Roman'))",sep="")
  step10 <- paste("lm_b1a_popgov_",i," <- lm(",i,"_demean ~ populism_y_demean, popopp_",i,")", sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5, step6, step7, step8, count, step9, step10)))
}

###Figure B1a
figureb1a <- grid.arrange(popopp_indlib_plot + ylab("Individual liberties") + xlab(""), popopp_ruleoflaw_plot + ylab("Rule of law") + xlab(""), popopp_public_plot + ylab("Public sphere") + xlab(""), popopp_compet_plot + ylab("Competition") + xlab(""), popopp_mutucons_plot + ylab("Mutual constraints") + xlab(""), popopp_govcap_plot + ylab("Governmental capacity") + xlab(""), popopp_transpar_plot + ylab("Transparency"), popopp_particip_plot + ylab("Participation"), popopp_repres_plot + ylab("Representation"), ncol=3,nrow=3)
ggsave(file="./figureb1a.jpg", figureb1a, width =  11.69, height = 8.27)

###Table B1a
stargazer(lm_b1a_popgov_indlib, lm_b1a_popgov_ruleoflaw, lm_b1a_popgov_public, lm_b1a_popgov_compet, lm_b1a_popgov_mutucons, lm_b1a_popgov_govcap, lm_b1a_popgov_transpar, lm_b1a_popgov_particip, lm_b1a_popgov_repres, title = "Table B1a: Populists in opposition and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years of populism", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")


##########Figure 3b / B1b: POPULISTS IN GOVERNMENT##########

for (i in (functions_var)) {
  step1 <- paste("popgov_",i," <- subset(af_wide, populism_gov == 1 & !is.na(",i,"))",sep="")
  step2 <- paste("popgov_",i," <- popgov_",i,"[c('gwid','gwid_populism_y','populism_y','",i,"')]",sep="")
  step3 <- paste("popgov_",i,"$",i,"_gm <- mean(popgov_",i,"$",i,", na.rm=T)",sep="")
  step4 <- paste("popgov_",i,"$populism_y_gm <- mean(popgov_",i,"$populism_y, na.rm=T)",sep="")
  step5 <- paste("popgov_",i," <- popgov_",i," %>% group_by(gwid_populism_y) %>% mutate(",i,"_periodm = mean(",i,",na.rm=T))",sep="")
  step6 <- paste("popgov_",i," <- popgov_",i," %>% group_by(gwid_populism_y) %>% mutate(populism_y_periodm = mean(populism_y,na.rm=T))",sep="")
  step7 <- paste("popgov_",i,"$",i,"_demean <- popgov_",i,"$",i," - popgov_",i,"$",i,"_periodm + popgov_",i,"$",i,"_gm",sep="")
  step8 <- paste("popgov_",i,"$populism_y_demean <- popgov_",i,"$populism_y - popgov_",i,"$populism_y_periodm + popgov_",i,"$populism_y_gm",sep="")
  count <- paste("popgov_",i,"$p10 <- ifelse(nrow(unique(popgov_",i,"[c('gwid')])) >= 10, 1, 0)", sep="")
  step9 <- paste("popgov_",i,"_plot <- ggplot() +
                 geom_line(data = popgov_",i,", aes(x = populism_y, y = ",i,", group = gwid_populism_y)) +
                 stat_smooth(data = popgov_",i,", aes(x = populism_y_demean, y = ",i,"_demean), method = 'lm', fullrange = T, xseq = seq(0, max(popgov_",i,"$p10 * popgov_",i,"$populism_y, na.rm=T), length=80), col = 'red') + 
                 labs(x ='Populism in gov. (y.)') +
                 theme(legend.position='none',plot.margin=unit(c(0.25,0.5,0.25,0.25), 'cm')) +
                 coord_cartesian(ylim = c(0, 100), xlim =c(1,10)) + 
                 scale_x_continuous(expand = c(0,0), breaks = c(2,4,6,8,10)) +
                 theme(text = element_text(size=rel(3.25))) +
                 theme(text=element_text(family='Times New Roman'))",sep="")
  step10 <- paste("lm_b1b_popgov_",i," <- lm(",i,"_demean ~ populism_y_demean, popgov_",i,")", sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5, step6, step7, step8, count, step9,step10)))
}

###Figure B1b
figureb1b <- grid.arrange(popgov_indlib_plot + ylab("Individual liberties") + xlab(""), popgov_ruleoflaw_plot + ylab("Rule of law") + xlab(""), popgov_public_plot + ylab("Public sphere") + xlab(""), popgov_compet_plot + ylab("Competition") + xlab(""), popgov_mutucons_plot + ylab("Mutual constraints") + xlab(""), popgov_govcap_plot + ylab("Governmental capacity") + xlab(""), popgov_transpar_plot + ylab("Transparency"), popgov_particip_plot + ylab("Participation"), popgov_repres_plot + ylab("Representation"), ncol=3,nrow=3)
ggsave(file="./figureb1b.jpg", figureb1b, width =  11.69, height = 8.27)

###Table B1b
stargazer(lm_b1b_popgov_indlib, lm_b1b_popgov_ruleoflaw, lm_b1b_popgov_public, lm_b1b_popgov_compet, lm_b1b_popgov_mutucons, lm_b1b_popgov_govcap, lm_b1b_popgov_transpar, lm_b1b_popgov_particip, lm_b1b_popgov_repres, title = "Table B1b: Populist government and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years of populism", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")


##########Figure 3c / B1c: POPULISTS IN DOMINANT GOVERNMENT##########

for (i in (functions_var)) {
  step1 <- paste("popgovd_",i," <- subset(af_wide, populism_gov_dom == 1 & !is.na(",i,"))",sep="")
  step2 <- paste("popgovd_",i," <- popgovd_",i,"[c('gwid','gwid_populism_govd_y','populism_govd_y','",i,"')]",sep="")
  step3 <- paste("popgovd_",i,"$",i,"_gm <- mean(popgovd_",i,"$",i,", na.rm=T)",sep="")
  step4 <- paste("popgovd_",i,"$populism_y_gm <- mean(popgovd_",i,"$populism_govd_y, na.rm=T)",sep="")
  step5 <- paste("popgovd_",i," <- popgovd_",i," %>% group_by(gwid_populism_govd_y) %>% mutate(",i,"_periodm = mean(",i,",na.rm=T))",sep="")
  step6 <- paste("popgovd_",i," <- popgovd_",i," %>% group_by(gwid_populism_govd_y) %>% mutate(populism_y_periodm = mean(populism_govd_y,na.rm=T))",sep="")
  step7 <- paste("popgovd_",i,"$",i,"_demean <- popgovd_",i,"$",i," - popgovd_",i,"$",i,"_periodm + popgovd_",i,"$",i,"_gm",sep="")
  step8 <- paste("popgovd_",i,"$populism_y_demean <- popgovd_",i,"$populism_govd_y - popgovd_",i,"$populism_y_periodm + popgovd_",i,"$populism_y_gm",sep="")
  count <- paste("popgovd_",i,"$p10 <- ifelse(nrow(unique(popgovd_",i,"[c('gwid')])) >= 10, 1, 0)", sep="")
  step9 <- paste("popgovd_",i,"_plot <- ggplot() +
                 geom_line(data = popgovd_",i,", aes(x = populism_govd_y, y = ",i,", group = gwid_populism_govd_y)) +
                 stat_smooth(data = popgovd_",i,", aes(x = populism_y_demean, y = ",i,"_demean), method = 'lm', fullrange=TRUE, , xseq = seq(0, max(popgovd_",i,"$p10 * popgovd_",i,"$populism_govd_y, na.rm=T), length=80), col = 'red') + 
                 labs(x ='Populism in gov. (dom., y.)') +
                 theme(legend.position='none',plot.margin=unit(c(0.25,0.5,0.25,0.25), 'cm')) +
                 coord_cartesian(ylim = c(0, 100), xlim =c(1,10)) + 
                 scale_x_continuous(expand = c(0,0), lim = c(0,10), breaks = c(2,4,6,8,10)) +
                 theme(text = element_text(size=rel(3.25))) +
                 theme(text=element_text(family='Times New Roman'))",sep="")
  step10 <- paste("lm_b1c_popgov_",i," <- lm(",i,"_demean ~ populism_y_demean, popgovd_",i,")", sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5, step6, step7, step8, count, step9, step10)))
}

###Figure B1c
figureb1c <- grid.arrange(popgovd_indlib_plot + ylab("Individual liberties") + xlab(""), popgovd_ruleoflaw_plot + ylab("Rule of law") + xlab(""), popgovd_public_plot + ylab("Public sphere") + xlab(""), popgovd_compet_plot + ylab("Competition") + xlab(""), popgovd_mutucons_plot + ylab("Mutual constraints") + xlab(""), popgovd_govcap_plot + ylab("govdernmental capacity") + xlab(""), popgovd_transpar_plot + ylab("Transparency"), popgovd_particip_plot + ylab("Participation"), popgovd_repres_plot + ylab("Representation"), ncol=3,nrow=3)
ggsave(file="./figureb1c.jpg", figureb1c, width =  11.69, height = 8.27)

###Table B1c
stargazer(lm_b1c_popgov_indlib, lm_b1c_popgov_ruleoflaw, lm_b1c_popgov_public, lm_b1c_popgov_compet, lm_b1c_popgov_mutucons, lm_b1c_popgov_govcap, lm_b1c_popgov_transpar, lm_b1c_popgov_particip, lm_b1c_popgov_repres, title = "Table B1c: Dominant populist parties and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years of populism", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")


##########FIGURE 4a / B2a: EU ASSOCIATION AGREEMENT##########

for (i in (functions_var)) {
  step1 <- paste("bs_aa_",i," <- subset(af_wide, !is.na(bs_aa_yP) & !is.na(",i,"))",sep="")
  step2 <- paste("bs_aa_",i," <- bs_aa_",i,"[c('gwid','bs_aa_yP','",i,"')]",sep="")
  step3 <- paste("bs_aa_",i,"$",i,"_gm <- mean(bs_aa_",i,"$",i,", na.rm=T)",sep="")
  step4 <- paste("bs_aa_",i,"$bs_aa_yP_gm <- mean(bs_aa_",i,"$bs_aa_yP, na.rm=T)",sep="")
  step5 <- paste("bs_aa_",i," <- bs_aa_",i," %>% group_by(gwid) %>% mutate(",i,"_periodm = mean(",i,",na.rm=T))",sep="")
  step6 <- paste("bs_aa_",i," <- bs_aa_",i," %>% group_by(gwid) %>% mutate(bs_aa_yP_periodm = mean(bs_aa_yP,na.rm=T))",sep="")
  step7 <- paste("bs_aa_",i,"$",i,"_demean <- bs_aa_",i,"$",i," - bs_aa_",i,"$",i,"_periodm + bs_aa_",i,"$",i,"_gm",sep="")
  step8 <- paste("bs_aa_",i,"$bs_aa_yP_demean <- bs_aa_",i,"$bs_aa_yP - bs_aa_",i,"$bs_aa_yP_periodm + bs_aa_",i,"$bs_aa_yP_gm",sep="")
  countbef <- paste("bs_aa_",i,"$before10c <- ifelse(nrow(unique(subset(bs_aa_",i,", bs_aa_yP <= 0)[c('gwid')])) >= 10, 1, 0)", sep="")
  countaft <- paste("bs_aa_",i,"$after10c <- ifelse(nrow(unique(subset(bs_aa_",i,", bs_aa_yP >= 0)[c('gwid')])) >= 10, 1, 0)", sep="")
  step9 <- paste("bs_aa_",i,"_plot <- ggplot() +
                 geom_line(data = subset(bs_aa_",i,", bs_aa_yP >= 0), aes(x = bs_aa_yP, y = ",i,", group = gwid, linetype = '1')) +
                 stat_smooth(data = subset(bs_aa_",i,", bs_aa_yP_demean >= 0), aes(x = bs_aa_yP_demean, y = ",i,"_demean, linetype = '1'), method = 'lm', fullrange=T, xseq = seq(0, max(bs_aa_",i,"$after10c * bs_aa_",i,"$bs_aa_yP), length=80), col = 'red') +
                 geom_line(data = subset(bs_aa_",i,", bs_aa_yP <= 0), aes(x = bs_aa_yP, y = ",i,", group = gwid, linetype = '2')) +
                 stat_smooth(data = subset(bs_aa_",i,", bs_aa_yP_demean <= 0), aes(x = bs_aa_yP_demean, y = ",i,"_demean, linetype = '2'), method = 'lm', fullrange=T, xseq = seq(min(bs_aa_",i,"$before10c * bs_aa_",i,"$bs_aa_yP), 0, length=80), col = 'red') + 
                 geom_vline(xintercept=0) +
                 labs(x ='Y. before/after\nassoc. agreement') +
                 theme(legend.position='none',plot.margin=unit(c(0.25,0.5,0.25,0.25), 'cm')) +
                 coord_cartesian(ylim = c(0, 100), xlim =c(-5,10)) + 
                 scale_x_continuous(expand = c(0,0), lim = c(-5,10), breaks = c(-4,-2,0,2,4,6,8,10)) +
                 theme(text = element_text(size=rel(4))) +
                 theme(text=element_text(family='Times New Roman'))",sep="")
  step10a <- paste("lm_b2aafter_popgov_",i," <- lm(",i,"_demean ~ bs_aa_yP_demean, subset = (bs_aa_yP_demean >= 0), bs_aa_",i,")", sep="")
  step10b <- paste("lm_b2abefore_popgov_",i," <- lm(",i,"_demean ~ bs_aa_yP_demean, subset = (bs_aa_yP_demean <= 0), bs_aa_",i,")", sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5, step6, step7, step8, countbef, countaft, step9, step10a, step10b)))
}
###changes
#labs(x ='Y. before/after\nassoc. agreement') +#\n -> two lines
#theme(text = element_text(size=rel(4))) + #4 instead of 3.25

###Figure B2a
figureb2a <- grid.arrange(bs_aa_indlib_plot + ylab("Individual liberties") + xlab(""), bs_aa_ruleoflaw_plot + ylab("Rule of law") + xlab(""), bs_aa_public_plot + ylab("Public sphere") + xlab(""), bs_aa_compet_plot + ylab("Competition") + xlab(""), bs_aa_mutucons_plot + ylab("Mutual constraints") + xlab(""), bs_aa_govcap_plot + ylab("Governmental capacity") + xlab(""), bs_aa_transpar_plot + ylab("Transparency"), bs_aa_particip_plot + ylab("Participation"), bs_aa_repres_plot + ylab("Representation"), ncol=3,nrow=3)
ggsave(file="./figureb2a.jpg", figureb2a, width =  11.69, height = 8.27)

###Tables B2a (1), B2a (2)
stargazer(lm_b2abefore_popgov_indlib, lm_b2abefore_popgov_ruleoflaw, lm_b2abefore_popgov_public, lm_b2abefore_popgov_compet, lm_b2abefore_popgov_mutucons, lm_b2abefore_popgov_govcap, lm_b2abefore_popgov_transpar, lm_b2abefore_popgov_particip, lm_b2abefore_popgov_repres, title = "Table B2a (1): Years before EU association agreements and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years before association agreement", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")
stargazer(lm_b2aafter_popgov_indlib, lm_b2aafter_popgov_ruleoflaw, lm_b2aafter_popgov_public, lm_b2aafter_popgov_compet, lm_b2aafter_popgov_mutucons, lm_b2aafter_popgov_govcap, lm_b2aafter_popgov_transpar, lm_b2aafter_popgov_particip, lm_b2aafter_popgov_repres, title = "Table B2a (2): Years after EU association agreements and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years after association agreement", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")


##########FIGURE 4b / B2b: EU WEAKLY CREDIBLE ACCESSION PERSPECTIVE##########

for (i in (functions_var)) {
  step1 <- paste("bs_wcap_",i," <- subset(af_wide, !is.na(bs_wcap_yP) & !is.na(",i,"))",sep="")
  step2 <- paste("bs_wcap_",i," <- bs_wcap_",i,"[c('gwid','bs_wcap_yP','",i,"')]",sep="")
  step3 <- paste("bs_wcap_",i,"$",i,"_gm <- mean(bs_wcap_",i,"$",i,", na.rm=T)",sep="")
  step4 <- paste("bs_wcap_",i,"$bs_wcap_yP_gm <- mean(bs_wcap_",i,"$bs_wcap_yP, na.rm=T)",sep="")
  step5 <- paste("bs_wcap_",i," <- bs_wcap_",i," %>% group_by(gwid) %>% mutate(",i,"_periodm = mean(",i,",na.rm=T))",sep="")
  step6 <- paste("bs_wcap_",i," <- bs_wcap_",i," %>% group_by(gwid) %>% mutate(bs_wcap_yP_periodm = mean(bs_wcap_yP,na.rm=T))",sep="")
  step7 <- paste("bs_wcap_",i,"$",i,"_demean <- bs_wcap_",i,"$",i," - bs_wcap_",i,"$",i,"_periodm + bs_wcap_",i,"$",i,"_gm",sep="")
  step8 <- paste("bs_wcap_",i,"$bs_wcap_yP_demean <- bs_wcap_",i,"$bs_wcap_yP - bs_wcap_",i,"$bs_wcap_yP_periodm + bs_wcap_",i,"$bs_wcap_yP_gm",sep="")
  step8.5 <- paste("bs_wcap_",i,"$beforeafter <- cut(bs_wcap_",i,"$bs_wcap_yP , breaks = c(-Inf, 0.00001, Inf), labels = c('before', 'after'))", sep="")
  countbef <- paste("bs_wcap_",i,"$before10c <- ifelse(nrow(unique(subset(bs_wcap_",i,", bs_wcap_yP <= 0)[c('gwid')])) >= 10, 1, 0)", sep="")
  countaft <- paste("bs_wcap_",i,"$after10c <- ifelse(nrow(unique(subset(bs_wcap_",i,", bs_wcap_yP >= 0)[c('gwid')])) >= 10, 1, 0)", sep="")
  step9 <- paste("bs_wcap_",i,"_plot <- ggplot() +
                 geom_line(data = subset(bs_wcap_",i,", bs_wcap_yP >= 0), aes(x = bs_wcap_yP, y = ",i,", group = gwid, linetype = '1')) +
                 stat_smooth(data = subset(bs_wcap_",i,", bs_wcap_yP_demean >= 0), aes(x = bs_wcap_yP_demean, y = ",i,"_demean, linetype = '1'), method = 'lm', fullrange=T, xseq = seq(0,max(bs_wcap_",i,"$after10c * bs_wcap_",i,"$bs_wcap_yP), length=80), col = 'red') + 
                 geom_line(data = subset(bs_wcap_",i,", bs_wcap_yP <= 0), aes(x = bs_wcap_yP, y = ",i,", group = gwid, linetype = '2')) +
                 stat_smooth(data = subset(bs_wcap_",i,", bs_wcap_yP_demean <= 0), aes(x = bs_wcap_yP_demean, y = ",i,"_demean, linetype = '2'), method = 'lm', fullrange=T, xseq = seq(min(bs_wcap_",i,"$after10c * bs_wcap_",i,"$bs_wcap_yP),0, length=80), col = 'red') +
                 geom_vline(xintercept=0) +
                 labs(x ='Y. before/after\nassoc. agreement') +
                 theme(legend.position='none',plot.margin=unit(c(0.25,0.5,0.25,0.25), 'cm')) +
                 coord_cartesian(ylim = c(0, 100), xlim =c(-5,10)) + 
                 scale_x_continuous(expand = c(0,0), lim = c(-5,10), breaks = c(-4,-2,0,2,4,6,8,10)) +
                 theme(text = element_text(size=rel(4))) +
                 theme(text=element_text(family='Times New Roman'))",sep="")
  step10a <- paste("lm_b2bafter_",i," <- lm(",i,"_demean ~ bs_wcap_yP_demean, subset = (bs_wcap_yP_demean >= 0), bs_wcap_",i,")", sep="")
  step10b <- paste("lm_b2bbefore_",i," <- lm(",i,"_demean ~ bs_wcap_yP_demean, subset = (bs_wcap_yP_demean <= 0), bs_wcap_",i,")", sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5, step6, step7, step8, countbef, countaft, step9, step10a, step10b)))
}

###Figure B2b
figureb2b <- grid.arrange(bs_wcap_indlib_plot + ylab("Individual liberties") + xlab(""), bs_wcap_ruleoflaw_plot + ylab("Rule of law") + xlab(""), bs_wcap_public_plot + ylab("Public sphere") + xlab(""), bs_wcap_compet_plot + ylab("Competition") + xlab(""), bs_wcap_mutucons_plot + ylab("Mutual constraints") + xlab(""), bs_wcap_govcap_plot + ylab("Governmental capacity") + xlab(""), bs_wcap_transpar_plot + ylab("Transparency"), bs_wcap_particip_plot + ylab("Participation"), bs_wcap_repres_plot + ylab("Representation"), ncol=3,nrow=3)
ggsave(file="./figureb2b.jpg", figureb2b, width =  11.69, height = 8.27)

###Tables B2b (1), B2b (2)
stargazer(lm_b2bbefore_indlib, lm_b2bbefore_ruleoflaw, lm_b2bbefore_public, lm_b2bbefore_compet, lm_b2bbefore_mutucons, lm_b2bbefore_govcap, lm_b2bbefore_transpar, lm_b2bbefore_particip, lm_b2bbefore_repres, title = "Table B2b (1): Years before weakly credible accession perspectives and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years before weakly credible accession perspective", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")
stargazer(lm_b2bafter_indlib, lm_b2bafter_ruleoflaw, lm_b2bafter_public, lm_b2bafter_compet, lm_b2bafter_mutucons, lm_b2bafter_govcap, lm_b2bafter_transpar, lm_b2bafter_particip, lm_b2bafter_repres, title = "Table B2b (2): Years after weakly credible accession perspectives and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years after weakly credible accession perspective", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")


##########Figure 4c / B2c: EU CREDIBLE ACCESSION PERSPECTIVE##########

for (i in (functions_var)) {
  step1 <- paste("bs_cap_",i," <- subset(af_wide, !is.na(bs_cap_yP) & !is.na(",i,"))",sep="")
  step2 <- paste("bs_cap_",i," <- bs_cap_",i,"[c('gwid','bs_cap_yP','",i,"')]",sep="")
  step3 <- paste("bs_cap_",i,"$",i,"_gm <- mean(bs_cap_",i,"$",i,", na.rm=T)",sep="")
  step4 <- paste("bs_cap_",i,"$bs_cap_yP_gm <- mean(bs_cap_",i,"$bs_cap_yP, na.rm=T)",sep="")
  step5 <- paste("bs_cap_",i," <- bs_cap_",i," %>% group_by(gwid) %>% mutate(",i,"_periodm = mean(",i,",na.rm=T))",sep="")
  step6 <- paste("bs_cap_",i," <- bs_cap_",i," %>% group_by(gwid) %>% mutate(bs_cap_yP_periodm = mean(bs_cap_yP,na.rm=T))",sep="")
  step7 <- paste("bs_cap_",i,"$",i,"_demean <- bs_cap_",i,"$",i," - bs_cap_",i,"$",i,"_periodm + bs_cap_",i,"$",i,"_gm",sep="")
  step8 <- paste("bs_cap_",i,"$bs_cap_yP_demean <- bs_cap_",i,"$bs_cap_yP - bs_cap_",i,"$bs_cap_yP_periodm + bs_cap_",i,"$bs_cap_yP_gm",sep="")
  step8.5 <- paste("bs_cap_",i,"$beforeafter <- cut(bs_cap_",i,"$bs_cap_yP , breaks = c(-Inf, 0.00001, Inf), labels = c('before', 'after'))", sep="")
  countbef <- paste("bs_cap_",i,"$before10c <- ifelse(nrow(unique(subset(bs_cap_",i,", bs_cap_yP <= 0)[c('gwid')])) >= 10, 1, 0)", sep="")
  countaft <- paste("bs_cap_",i,"$after10c <- ifelse(nrow(unique(subset(bs_cap_",i,", bs_cap_yP >= 0)[c('gwid')])) >= 10, 1, 0)", sep="")
  step9 <- paste("bs_cap_",i,"_plot <- ggplot() +
                 geom_line(data = subset(bs_cap_",i,", bs_cap_yP >= 0), aes(x = bs_cap_yP, y = ",i,", group = gwid, linetype = '1')) +
                 stat_smooth(data = subset(bs_cap_",i,", bs_cap_yP_demean >= 0), aes(x = bs_cap_yP_demean, y = ",i,"_demean, linetype = '1'), method = 'lm', fullrange=T, xseq = seq(0,max(bs_cap_",i,"$after10c * bs_cap_",i,"$bs_cap_yP), length=80), col = 'red') + 
                 geom_line(data = subset(bs_cap_",i,", bs_cap_yP <= 0), aes(x = bs_cap_yP, y = ",i,", group = gwid, linetype = '2')) +
                 stat_smooth(data = subset(bs_cap_",i,", bs_cap_yP_demean <= 0), aes(x = bs_cap_yP_demean, y = ",i,"_demean, linetype = '2'), method = 'lm', fullrange=T, xseq = seq(min(bs_cap_",i,"$after10c * bs_cap_",i,"$bs_cap_yP),0, length=80), col = 'red') + 
                 geom_vline(xintercept=0) +
                 labs(x ='Y. before/after\ncredible acc.') +
                 theme(legend.position='none',plot.margin=unit(c(0.25,0.5,0.25,0.25), 'cm')) +
                 coord_cartesian(ylim = c(0, 100), xlim =c(-5,10)) + 
                 scale_x_continuous(expand = c(0,0), lim = c(-5,10), breaks = c(-4,-2,0,2,4,6,8,10)) +
                 theme(text = element_text(size=rel(4))) +
                 theme(text=element_text(family='Times New Roman'))",sep="")
  step10a <- paste("lm_b2cafter_",i," <- lm(",i,"_demean ~ bs_cap_yP, subset = (bs_cap_yP >= 0), bs_cap_",i,")", sep="")
  step10b <- paste("lm_b2cbefore_",i," <- lm(",i,"_demean ~ bs_cap_yP, subset = (bs_cap_yP <= 0), bs_cap_",i,")", sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5, step6, step7, step8, countbef, countaft, step9, step10a, step10b)))
}

###Figure B2c
figureb2c <- grid.arrange(bs_cap_indlib_plot + ylab("Individual liberties") + xlab(""), bs_cap_ruleoflaw_plot + ylab("Rule of law") + xlab(""), bs_cap_public_plot + ylab("Public sphere") + xlab(""), bs_cap_compet_plot + ylab("Competition") + xlab(""), bs_cap_mutucons_plot + ylab("Mutual constraints") + xlab(""), bs_cap_govcap_plot + ylab("Governmental capacity") + xlab(""), bs_cap_transpar_plot + ylab("Transparency"), bs_cap_particip_plot + ylab("Participation"), bs_cap_repres_plot + ylab("Representation"), ncol=3,nrow=3)
ggsave(file="./figureb2c.jpg", figureb2c, width =  11.69, height = 8.27)

###Tables B2c (1), B2c (2)
stargazer(lm_b2cbefore_indlib, lm_b2cbefore_ruleoflaw, lm_b2cbefore_public, lm_b2cbefore_compet, lm_b2cbefore_mutucons, lm_b2cbefore_govcap, lm_b2cbefore_transpar, lm_b2cbefore_particip, lm_b2cbefore_repres, title = "Table B2c (1): Years before credible accession perspectives and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years before credible accession perspective", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")
stargazer(lm_b2cafter_indlib, lm_b2cafter_ruleoflaw, lm_b2cafter_public, lm_b2cafter_compet, lm_b2cafter_mutucons, lm_b2cafter_govcap, lm_b2cafter_transpar, lm_b2cafter_particip, lm_b2cafter_repres, title = "Table B2c (2): Years after credible accession perspectives and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years after credible accession perspective", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")


##########FIGURE 4D / B2d: EU MEMBERSHIP##########

for (i in (functions_var)) {
  step1 <- paste("bs_eum_",i," <- subset(af_wide, !is.na(bs_eum_yP) & !is.na(",i,"))",sep="")
  step2 <- paste("bs_eum_",i," <- bs_eum_",i,"[c('gwid','bs_eum_yP','",i,"')]",sep="")
  step3 <- paste("bs_eum_",i,"$",i,"_gm <- mean(bs_eum_",i,"$",i,", na.rm=T)",sep="")
  step4 <- paste("bs_eum_",i,"$bs_eum_yP_gm <- mean(bs_eum_",i,"$bs_eum_yP, na.rm=T)",sep="")
  step5 <- paste("bs_eum_",i," <- bs_eum_",i," %>% group_by(gwid) %>% mutate(",i,"_periodm = mean(",i,",na.rm=T))",sep="")
  step6 <- paste("bs_eum_",i," <- bs_eum_",i," %>% group_by(gwid) %>% mutate(bs_eum_yP_periodm = mean(bs_eum_yP,na.rm=T))",sep="")
  step7 <- paste("bs_eum_",i,"$",i,"_demean <- bs_eum_",i,"$",i," - bs_eum_",i,"$",i,"_periodm + bs_eum_",i,"$",i,"_gm",sep="")
  step8 <- paste("bs_eum_",i,"$bs_eum_yP_demean <- bs_eum_",i,"$bs_eum_yP - bs_eum_",i,"$bs_eum_yP_periodm + bs_eum_",i,"$bs_eum_yP_gm",sep="")
  step8.5 <- paste("bs_eum_",i,"$beforeafter <- cut(bs_eum_",i,"$bs_eum_yP , breaks = c(-Inf, 0.00001, Inf), labels = c('before', 'after'))", sep="")
  countbef <- paste("bs_eum_",i,"$before10c <- ifelse(nrow(unique(subset(bs_eum_",i,", bs_eum_yP <= 0)[c('gwid')])) >= 10, 1, 0)", sep="")
  countaft <- paste("bs_eum_",i,"$after10c <- ifelse(nrow(unique(subset(bs_eum_",i,", bs_eum_yP >= 0)[c('gwid')])) >= 10, 1, 0)", sep="")
  step9 <- paste("bs_eum_",i,"_plot <- ggplot() +
                 geom_line(data = subset(bs_eum_",i,", bs_eum_yP >= 0), aes(x = bs_eum_yP, y = ",i,", group = gwid, linetype = '1')) +
                 stat_smooth(data = subset(bs_eum_",i,", bs_eum_yP_demean >= 0), aes(x = bs_eum_yP_demean, y = ",i,"_demean, linetype = '1'), method = 'lm', fullrange=T, xseq = seq(0,max(bs_eum_",i,"$after10c * bs_eum_",i,"$bs_eum_yP), length=80), col = 'red') + 
                 geom_line(data = subset(bs_eum_",i,", bs_eum_yP <= 0), aes(x = bs_eum_yP, y = ",i,", group = gwid, linetype = '2')) +
                 stat_smooth(data = subset(bs_eum_",i,", bs_eum_yP_demean <= 0), aes(x = bs_eum_yP_demean, y = ",i,"_demean, linetype = '2'), method = 'lm', fullrange=T, xseq = seq(min(bs_eum_",i,"$after10c * bs_eum_",i,"$bs_eum_yP),0, length=80), col = 'red') + 
                 geom_vline(xintercept=0) +
                 labs(x ='Y. before/after\nmembership') +
                 theme(legend.position='none',plot.margin=unit(c(0.25,0.5,0.25,0.25), 'cm')) +
                 coord_cartesian(ylim = c(0, 100), xlim =c(-5,10)) + 
                 scale_x_continuous(expand = c(0,0), lim = c(-5,10), breaks = c(-4,-2,0,2,4,6,8,10)) +
                 theme(text = element_text(size=rel(4))) +
                 theme(text=element_text(family='Times New Roman'))",sep="")
  step10a <- paste("lm_b2dafter_",i," <- lm(",i,"_demean ~ bs_eum_yP_demean, subset = (bs_eum_yP_demean >= 0), bs_eum_",i,")", sep="")
  step10b <- paste("lm_b2dbefore_",i," <- lm(",i,"_demean ~ bs_eum_yP_demean, subset = (bs_eum_yP_demean <= 0), bs_eum_",i,")", sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5, step6, step7, step8, countbef, countaft, step9, step10a, step10b)))
}

###Figure B2d
figureb2d <- grid.arrange(bs_eum_indlib_plot + ylab("Individual liberties") + xlab(""), bs_eum_ruleoflaw_plot + ylab("Rule of law") + xlab(""), bs_eum_public_plot + ylab("Public sphere") + xlab(""), bs_eum_compet_plot + ylab("Competition") + xlab(""), bs_eum_mutucons_plot + ylab("Mutual constraints") + xlab(""), bs_eum_govcap_plot + ylab("Governmental capacity") + xlab(""), bs_eum_transpar_plot + ylab("Transparency"), bs_eum_particip_plot + ylab("Participation"), bs_eum_repres_plot + ylab("Representation"), ncol=3,nrow=3)
ggsave(file="./figureb2d.jpg", figureb2d, width =  11.69, height = 8.27)

###Tables B2d (1), B2d (2)
stargazer(lm_b2dbefore_indlib, lm_b2dbefore_ruleoflaw, lm_b2dbefore_public, lm_b2dbefore_compet, lm_b2dbefore_mutucons, lm_b2dbefore_govcap, lm_b2dbefore_transpar, lm_b2dbefore_particip, lm_b2dbefore_repres, title = "Table B2d (1): Years before EU accession  and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years before EU accession", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")
stargazer(lm_b2dafter_indlib, lm_b2dafter_ruleoflaw, lm_b2dafter_public, lm_b2dafter_compet, lm_b2dafter_mutucons, lm_b2dafter_govcap, lm_b2dafter_transpar, lm_b2dafter_particip, lm_b2dafter_repres, title = "Table B2d (2): Years after EU accession and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years after EU accession", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")


##########FIGURE 5 / B3: FINANCIAL CRISIS##########
for (i in (functions_var)) {
  step1 <- paste("crisis_",i," <- subset(af_wide, !is.na(",i,") & year >= 2003)",sep="")
  step2 <- paste("crisis_",i," <- crisis_",i,"[c('gwid','year','",i,"')]",sep="")
  step3 <- paste("crisis_",i,"$",i,"_gm <- mean(crisis_",i,"$",i,", na.rm=T)",sep="")
  step4 <- paste("crisis_",i,"$year_gm <- mean(crisis_",i,"$year, na.rm=T)",sep="")
  step5 <- paste("crisis_",i," <- crisis_",i," %>% group_by(gwid) %>% mutate(",i,"_periodm = mean(",i,",na.rm=T))",sep="")
  step6 <- paste("crisis_",i," <- crisis_",i," %>% group_by(gwid) %>% mutate(year_periodm = mean(year,na.rm=T))",sep="")
  step7 <- paste("crisis_",i,"$",i,"_demean <- crisis_",i,"$",i," - crisis_",i,"$",i,"_periodm + crisis_",i,"$",i,"_gm",sep="")
  step8 <- paste("crisis_",i,"$year_demean <- crisis_",i,"$year - crisis_",i,"$year_periodm + crisis_",i,"$year_gm",sep="")
  step9 <- paste("crisis_",i,"_plot <- ggplot() +
                 geom_line(data = subset(crisis_",i,", year >= 2008), aes(x = year, y = ",i,", group = gwid, linetype = '1')) +
                 stat_smooth(data = subset(crisis_",i,", year_demean >= 2008), aes(x = year_demean, y = ",i,"_demean, linetype = '1'), method = 'lm', fullrange=T, xseq = seq(2008,2013, length=80), col = 'red') + 
                 geom_line(data = subset(crisis_",i,", year <= 2008), aes(x = year, y = ",i,", group = gwid, linetype = '2')) +
                 stat_smooth(data = subset(crisis_",i,", 2003 & year_demean <= 2008), aes(x = year_demean, y = ",i,"_demean, linetype = '2'), method = 'lm', fullrange=T, xseq = seq(2003,2008, length=80), col = 'red') + 
                 geom_vline(xintercept=2008) +
                 labs(x ='Y. before/after financial crisis') +
                 theme(legend.position='none',plot.margin=unit(c(0.25,0.5,0.25,0.25), 'cm')) +
                 coord_cartesian(ylim = c(0, 100), xlim =c(2003,2013)) + 
                 scale_x_continuous(expand = c(0,0), lim = c(2003,2013), breaks = c(2003,2005,2007,2009,2011,2013)) +
                 theme(text = element_text(size=rel(3.25))) +
                 theme(text=element_text(family='Times New Roman'))",sep="")
  step10a <- paste("lm_b3after_",i," <- lm(",i,"_demean ~ year_demean, subset = (year_demean >= 2008), crisis_",i,")", sep="")
  step10b <- paste("lm_b3before_",i," <- lm(",i,"_demean ~ year_demean, subset = (year_demean <= 2008), crisis_",i,")", sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5, step6, step7, step8, step9, step10a, step10b)))
}

###Figure B3
figureb3 <- grid.arrange(crisis_indlib_plot + ylab("Individual liberties") + xlab(""), crisis_ruleoflaw_plot + ylab("Rule of law") + xlab(""), crisis_public_plot + ylab("Public sphere") + xlab(""), crisis_compet_plot + ylab("Competition") + xlab(""), crisis_mutucons_plot + ylab("Mutual constraints") + xlab(""), crisis_govcap_plot + ylab("Governmental capacity") + xlab(""), crisis_transpar_plot + ylab("Transparency"), crisis_particip_plot + ylab("Participation"), crisis_repres_plot + ylab("Representation"), ncol=3,nrow=3)
ggsave(file="./figureb3.jpg", figureb3, width =  11.69, height = 8.27)

###Tables B3 (1), B3 (2)
stargazer(lm_b3before_indlib, lm_b3before_ruleoflaw, lm_b3before_public, lm_b3before_compet, lm_b3before_mutucons, lm_b3before_govcap, lm_b3before_transpar, lm_b3before_particip, lm_b3before_repres, title = "Table B3 (1): Years before financial crisis and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years before financial crisis", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")
stargazer(lm_b3after_indlib, lm_b3after_ruleoflaw, lm_b3after_public, lm_b3after_compet, lm_b3after_mutucons, lm_b3after_govcap, lm_b3after_transpar, lm_b3after_particip, lm_b3after_repres, title = "Table B3 (2): Years after financial crisis and trends in the quality of democracy (models used for figures).", colnames = F, dep.var.labels =c("Individual liberties", "Rule of law", "Public sphere", "Competitition", "Mutual constraints", "Governmental capacity", "Transparency", "Participation", "Representation"), covariate.labels = c("Years after financial crisis", "Constant"), type = "text",  notes = "All models use demeaned independent and dependent variables, analogously to a manual (populism-period) fixed-effect model.")


##########Main figures for article##########

figure3 <- grid.arrange(popopp_indlib_plot + ylab ("Individual liberties") + xlab(""), popgov_indlib_plot + ylab ("") + xlab(""), popgovd_indlib_plot + ylab ("") + xlab(""), popopp_ruleoflaw_plot + ylab ("Rule of law") + xlab(""), popgov_ruleoflaw_plot + ylab ("") + xlab(""), popgovd_ruleoflaw_plot + ylab ("") + xlab(""), popopp_compet_plot + ylab ("Competition") + xlab(""), popgov_compet_plot + ylab ("") + xlab(""), popgovd_compet_plot + ylab ("") + xlab(""), popopp_transpar_plot + ylab ("Transparency"), popgov_transpar_plot + ylab (""), popgovd_transpar_plot + ylab (""), ncol=3,nrow=4)
ggsave(file="./figure3.png", figure3, width =  8.27, height = 11.69, dpi = 600)
figure4 <- grid.arrange(bs_aa_indlib_plot + ylab ("Ind. liberties") + xlab(""), bs_wcap_indlib_plot + ylab ("") + xlab(""), bs_cap_indlib_plot + ylab ("") + xlab(""), bs_eum_indlib_plot + ylab ("") + xlab(""), bs_aa_ruleoflaw_plot + ylab ("Rule of law") + xlab(""), bs_wcap_ruleoflaw_plot + ylab ("") + xlab(""), bs_cap_ruleoflaw_plot + ylab ("") + xlab(""), bs_eum_ruleoflaw_plot + ylab ("") + xlab(""), bs_aa_compet_plot + ylab ("Competition") + xlab(""), bs_wcap_compet_plot + ylab ("") + xlab(""), bs_cap_compet_plot + ylab ("") + xlab(""), bs_eum_compet_plot + ylab ("") + xlab(""), bs_aa_transpar_plot + ylab ("Transparency"), bs_wcap_transpar_plot + ylab (""), bs_cap_transpar_plot + ylab (""), bs_eum_transpar_plot + ylab (""), ncol=4,nrow=4)
ggsave(file="./figure4.png", figure4, width =  11.69, height = 8.27, dpi = 600)
figure5 <- grid.arrange(crisis_indlib_plot + ylab("Individual liberties") + xlab("Y. before/after financial crisis"), crisis_ruleoflaw_plot + ylab("Rule of law") + xlab("Y. before/after financial crisis"), crisis_compet_plot + ylab("Competition") + xlab("Y. before/after financial crisis"), crisis_transpar_plot + ylab("Transparency") + xlab("Y. before/after financial crisis"), ncol=4,nrow=1)
ggsave(file="./figure5.png", figure5, width =  8.27, height = 2, dpi = 600)